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1 Introduction 

pydelay is a program which translates a system of delay differential equations (DDEs) into simulation C-code and 
compiles and runs the code (using scipy weave). This way it is easy to quickly implement a system of DDEs but 
you still have the speed of C. The Homepage can be found here: 

http://pydelay.sourceforge.net/ 

It is largely inspired by PyDSTool. 

The algorithm used is based on the Bogacki-Shampine method ' which is also implemented in Matlab's dde23 2 . 

We also want to mention PyDDE - a different python program for solving DDEs. 

License 

pydelay is licensed under the MIT License. 

1.1 Installation and requirements 

Unix: 

You need python, numpy and scipy and the gcc-compiler. To plot the solutions and run the examples you also 
need matplotlib. 

To install pydelay grab the latest tar.gz from the website and install the package in the usual way: 

cd pydelay-$version 
python setup. py install 

When the package is installed, you can get some info about the functions and the usage with: 

pydoc pydelay 

Windows: 

The solver has not been tested on a windows machine. It could perhaps work under cygwin. 

1 Bogacki, P. and Shampine, L. R, A 3(2) pair of Runge - Kutta formulas, Applied Mathematics Letters 2, 4, 321 ISSN 0893-9659, (1989). 

2 Shampine, L. F. and Thompson, S., Solving DDEs in Matlab, Appl. Num. Math. 37, 4, 441 ( 2001) 



1.2 An example 



The following example shows the basic usage. It solves the Mackey-Glass equations 3 for initial conditions which 
lead to a periodic orbit (see 4 for this example). 

# import pydelay and numpy and pylab 
import numpy as np 

import pylab as pi 

from pydelay import dde23 

# define the equations 
eqns = { 

'x' : '0.25 * x(t-tau) / (1.0 + pow (x (t-tau) , p) ) -0.1*x' 
} 

idefine the parameters 
params = { 

' tau' : 15, 

'p' : 10 

} 

# Initialise the solver 

dde = dde23 (eqns=eqns, params=params ) 

# set the simulation parameters 

# (solve from t=0 to t=1000 and limit the maximum step size to 1.0) 
dde . set_sim_params (tf inal=1000, dtmax=l . 0) 

# set the history of to the constant function 0.5 (using a python lambda function) 
histfunc = { 

' x' : lambda t : 0.5 
} 

dde . hist_f rom_f uncs (histfunc, 51) 

# run the simulator 
dde . run ( ) 

# Make a plot of x(t) vs x(t-tau): 

# Sample the solution twice with a stepsize of dt=0 . 1 : 

# once in the interval [515, 1000] 
soil = dde . sample (515, 1000, 0.1) 
xl = soil [' x' ] 

# and once between [500, 1000-15] 
sol2 = dde . sample (500, 1000-15, 0.1) 
x2 = sol2 [ ' x' ] 

pi . plot (xl , x2 ) 

pi .xlabel (' $x(t) $' ) 

pi . ylabel (' $x (t - 15) $' ) 

pi . show ( ) 



3 Mackey, M. C. and Glass, L. (1977). Pathological physiological conditions resulting from instabilities in physiological control system. 
Science, 197(4300):287-289. 

4 http://www.scholarpedia.org/article/Mackey-Glass_equation 
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2 Usage 

2.1 Defining the equations, delays and parameters 

Equations are defined using a python dictionary. The keys are the variable names and the entry is the right 
hand side of the differential equation. The string defining the equation has to be a valid C expression, i.e., use 
pow(a,b) instead of a* *b etc. 

Delays are written as (t -de lay) , where delay can be some expression involving parameters and numbers but 
not (yet) involving the time t or the dynamic variables: 

eqns = { 

'yl': '- yl * y2 (t-tau) + y2(t-1.0)', 
'y2': 'a * yl * y2 (t-2*tau) - y2' , 
'y3' : 'y2 - y2 (t- (tau+1) ) ' 

} 

Complex variables can be defined by adding ' : c ' or ' : C ' in the eqn-dictionary. The imaginary unit can be used 
through ' ii' in the equations: 

eqns = { 

'z:c': '(la + ii*wO + g*pow (abs ( z ) , 2 ) ) * z + b* ( z (t-tau) - z (t ) ) ' , 

} 

Parameters are defined in a separate dictionary where the keys are the parameter names, i.e.,: 



params = { 

'a' : 0.2, 



' tau' : 1.0 

} 



2.2 Setting the history 

The history of the variables is stored in the dictionary dde2 3 .hist. The keys are the variable names and there 
is an additional key ' t ' for the time array of the history. 

There is a second dictionary dde 2 3 . Vh i s t where the time derivatives of the history is stored (this is needed for 
the solver). When the solver is initialized, i.e.,: 

dde = dde23(eqns, params) 

the history of all variables (defined in eqns) is initialized to an array of length nn=101 filled with zeros. The 
time array is evenly spaced in the interval [-maxdelay, 0]. 

It is possible to manipulate these arrays directly, however this is not recommended since one easily ends up with 
an ill-defined history resulting for example in segfaults or false results. 

Instead use the following methods to set the history. 

hist_f rom_funcs (die, nn=101) 

Initialise the histories with the functions stored in the dictionary die. The keys are the variable names. The 
function will be called as f (t ) for t in [ -maxdelay, ] on nn samples in the interval. 

This function provides the simplest way to set the history. It is often convenient to use python lambda 
functions for f . This way you can define the history function in place. 

If any variable names are missing in the dictionaries, the history of these variables is set to zero and a 
warning is printed. If the dictionary contains keys not matching any variables these entries are ignored and 
a warning is printed. 

Example: Initialise the history of the variables x and y with cos and sin functions using a finer sampling 
resolution: 

from math import sin, cos 

histdic = { 

'x': lambda t: cos(0.2*t), 
'y': lambda t: sin(0.2*t) 

} 

dde . hist_f rom_f uncs (histdic, 500) 

hist_f rom_arrays (die, useend=True) 

Initialise the history using a dictionary of arrays with variable names as keys. Additionally a time array can 
be given corresponding to the key t. All arrays in die have to have the same lengths. 

If an array for t is given the history is interpreted as points (t , var ) . Otherwise the arrays will be evenly 
spaced out over the interval [-maxdelay, 0]. 

If useend is True the time array is shifted such that the end time is zero. This is useful if you want to use the 
result of a previous simulation as the history. 

If any variable names are missing in the dictionaries, the history of these variables is set to zero and a 
warning is printed. If the dictionary contains keys not matching any variables (or ' t ' ) these entries are 
ignored and a warning is printed. 

Example:: 



t = numpy . linspace ( , 1, 500) 
x = numpy . cos ( . 2 *t ) 
y = numpy . sin ( . 2 *t ) 



histdic 


= { 


' t' : 


t, 


' x' : 


x, 


,y, . 


y 


) 





dde . hist_f rom_arrays (histdic) 



Note that the previously used methods hist_f rom_dict, hist_f rom_array and hist_f rom_f unc (the 
last two without s) have been removed, since it was too easy to make mistakes with them. 



2.3 The solution 

After the solver has run, the solution (including the history) is stored in the dictionary dde 2 3 .sol. The keys are 
again the variable names and the time ' t ' . Since the solver uses an adaptive step size method, the solution is not 
sampled at regular times. 

To sample the solutions at regular (or other custom spaced) times there are two functions. 

sample (tstart=None, tfinal=None, dt=None) 

Sample the solution with dt steps between tstart and tfinal. 

tstart, tfinal Start and end value of the interval to sample. If nothing is specified tstart is set to zero and 
tfinal is set to the simulation end time. 

dt Sampling size used. If nothing is specified a reasonable value is calculated. 

Returns a dictionary with the sampled arrays. The keys are the variable names. The key ' t ' corresponds 
to the sampling times. 

sol_spl (t) 

Sample the solutions at times t. 

t Array of time points on which to sample the solution. 

Returns a dictionary with the sampled arrays. The keys are the variable names. The key ' t ' corresponds 
to the sampling times. 

These functions use a cubic spline interpolation of the solution data. 



2.4 Noise 

Noise can be included in the simulations. Note however, that the method used is quite crude (an Euler method 
will be added which is better suited for noise dominated dynamics). The deterministic terms are calculated with 
the usual Runge-Kutta method and then the noise term is added with the proper scaling of $ \ sqrt { dt } $ at the 
final step. To get accurate results one should use small time steps, i.e., dtmax should be set small enough. 

The noise is defined in a separate dictionary. The function gwn ( ) can be accessed in the noise string and is a 
Gaussian white noise term of unit variance. The following code specifies an Ornstein-Uhlenbeck process.: 

eqns = { ' x' : ' -x' } 

noise = { ' x' : ' D * gwn ( ) ' } 

params = { ' D' : 0.00001 } 

dde = dde23 (eqns=eqns, params=params , noise=noise) 

You can also use noise terms of other forms by specifying an appropriate C-function (see the section on custom 
C-code). 



2.5 Custom C-code 



You can access custom C-functions in your equations by adding the definition as supportcode for the solver. 
In the following example a function f ( w, t ) is defined through C-code and accessed in the eqn string.: 

# define the eqn f is the C- function defined below 
eqns = { 'x': '- x + k*x(t-tau) + A*f(w,t)' } 
params = { 

'k' : 0.1, 
'w' : 2.0, 
'A' : 0.5, 
'tau' : 10.0 

} 

mycode = " " " 

double f (double t, double w) { 
return sin (w * t ) ; 

} 

it M ii 

dde = dde23 (eqns=eqns, params=params , supportcode=mycode ) 

When defining custom code you have to be careful with the types. The type of complex variables in the C-code is 
Complex. Note in the above example that w has to be given as an input to the function, because the parameters 
can only be accessed from the eqns string and not inside the supportcode. (Should this be changed?) 

Using custom C-code is often useful for switching terms on and off. For example the Heaviside function may be 
defined and used as follows.: 

# define the eqn f is the C- function defined below 

eqns = { 'z:c': '(la+ii*w)*z - Heavi (t-t ) * K* (z-z (t-tau) ) ' } 
params = { 



K' 


. 1 


w' 


1.0, 


la' 


0.1, 


tau' 


Pi/ 


t0' 


2 *pi 



} 

mycode = " " " 

double Heavi (double t) { 
if (t>=0) 

return 1.0; 

else 

return 0.0; 

} 

ii ii ii 

dde = dde23 (eqns=eqns, params=params , supportcode=mycode ) 

This code would switch a control term on when t>t0. Note that Heavi (t-tO) does not get translated to a 
delay term, because Heavi is not a system variable. 

Since this scenario occurs so frequent the Heaviside function (as defined above) is included by default in the 
source code. 

2.6 Use and modify generated code 

The compilation of the generated code is done with scipy . weave. Instead of using weave to run the code you 
can directly access the generated code via the function dde2 3 . output_ccode ( ) . This function returns the 
generated code as a string which you can then store in a source file. 



To run the generated code manually you have to set the precompiler flag\ #def ine MANUAL (uncomment the 
line in the source file) to exclude the python related parts and include some other parts making the code a valid 
stand alone source file. After this the code should compile with g++ -lm -o prog source . cpp and you 
can run the program manually. 

You can specify the history of all variables in the source file by setting the for loops after the comment\ / * set 

the history here . . . */. 

Running the code manually can help you debug, if some problem occurs and also allows you to extend the code 
easily. 

2.7 Another example 

The following example shows some of the things discussed above. The code simulates the Lang-Kobayashi laser 
equations 5 

E'{t) = i(l + ia)nE + KE(t - r) 
Tn'(t)=p-n-(l + n)\E\ 2 

import numpy as np 
import pylab as pi 
from pydelay import dde23 

tfinal = 10000 
tau = 1000 

#the laser equations 
eqns = { 

'E:c': ' . 5* (1 . 0+ii*a) *E*n + K*E (t-tau) ' , 

'n' : ' (p - n - (1.0 +n) * pow (abs (E ) , 2 ) ) /T' 

} 



params = 




' a' 


4.0, 


'P' 


1.0, 


f <j> f 


200.0, 


' K' 


0.1, 


' tau' 


tau, 


' nu' 


10**-5, 


' n0' 


10 . 



} 

noise = { 'E': 'sqrt(0.5*nu*(n+n0)) * (gwn ( ) + ii*gwn ( ) ) ' } 

dde = dde23 (eqns=eqns, params=params , noise=noise) 
dde . set_sim_params (tf inal=tf inal) 

# use a dictionary to set the history 
thist = np . linspace ( , tau, tfinal) 
Ehist = np . zeros (len (thist) ) +1 . 
nhist = np . zeros (len (thist) ) -0 . 2 

die = {'t' : thist, ' E ' : Ehist, 'n': nhist} 

# 'useend' is True by default in hist_from_dict and thus the 

# time array is shifted correctly 
dde . hist_f rom_arrays (die) 

dde . run ( ) 



Lang, R. and Kobayashi, K. , External optical feedback effects on semiconductor injection laser properties, IEEE J. Quantum Electron. 
16, 347 (1980) 



t = dde.sol ['t' ] 
E = dde . sol [ ' E ' ] 
n = dde . sol [ ' n' ] 

spl = dde . sample (-tau, tfinal, 0.1) 

pi .plot (t [ : -1] , t[l:] - t[:-l], '0.8', label = ' step size') 

pi . plot ( spl [ ' t ' ] , abs ( spl [ ' E ' ] ) , ' g' , label = ' sampled solution') 

pi .plot (t, abs (E) , ' . ' , label = ' calculated points' ) 

pi . legend ( ) 

pi .xlabel (' $tS' ) 
pl.ylabel('$|E|$') 

pi . xlim ( (0 . 95*tf inal, tfinal)) 
pl.ylim( (0,3) ) 
pi . show ( ) 



1.5 




00 9600 9700 9800 9900 10000 

t 



3 Module Reference 



_init (eqns, params=None, noise=None, supportcode=", debug=False) 

Initialise the solver. 

eqns Dictionary defining for each variable the derivative. Delays are written as as (t 



example: 



eqns = { 

'yl': '- yl * y2(t-taul) + y2(t-tau2)', 
'y2': 'a * yl * y2(t-taul) - y2', 



' y3' : ' y2 - y2 (t-tau2) ' 
} 

You can also directly use numbers or combination of parameters as delays: 

eqns = { 

' xl' : ' -a*xl + xl (t - 1.0)', 
'x2': ' x2-b*xl (t-2 . 0*a+b) 
} 

At the moment only constant delays are supported. 

The string defining the equation has to be a valid C expression, i.e., use pow (a, b) instead of a**b 
etc. (this might change in the future): 

eqns = {'y': ' -2 . * sin(t) * pow (y (t-tau) , 2)'} 

Complex variable can be defined using : C or : c in the variable name. The imaginary unit can be used 
through ii in the equations: 

eqns = {'z:C: '(-la + ii * wO ) * z' } 
params Dictionary defining the parameters (including delays) used in eqns. example: 

params = { 

'a' : 1.0, 
' taul' : 1.0, 
'tau2' : 10.0 
} 

noise Dictionary for noise terms. The function gwn ( ) can be accessed in the noise string and provides a 
Gaussian white noise term of unit variance, example: 

noise = {'x': ' . 1 *gwn ( ) ' } 

debug If set to True the solver gives verbose output to stdout while running. 

set_sim_params (tfinal=100, AbsTol=9. 999999999999999 5 e-07, RelTol=0.001, 

dtmin=9.9999999999999995e-07, dtmax=None, dtO=None, Maxlter=1000000000.0) 

tfinal End time of the simulation (the simulation always starts at t=0). 

AbsTol, RelTol The relative and absolute error tolerance. If the estimated error e for a variable y obeys e 
<= AbsTol + RelTol*\y\ then the step is accepted. Otherwise the step will be repeated with a smaller 
step size. 

dtmin, dtmax Minimum and maximum step size used. 
dtO initial step size 

Maxlter maximum number of steps. The simulation stops if this is reached. 

hist_f rom_arrays (die, useend=True) 

Initialise the history using a dictionary of arrays with variable names as keys. Additionally a time array can 
be given corresponding to the key t. All arrays in die have to have the same lengths. 

If an array for t is given the history is interpreted as points (t , var ) . Otherwise the arrays will be evenly 
spaced out over the interval [-maxdelay, 0]. 

If useend is True the time array is shifted such that the end time is zero. This is useful if you want to use the 
result of a previous simulation as the history. 

If any variable names are missing in the dictionaries, the history of these variables is set to zero and a 
warning is printed. If the dictionary contains keys not matching any variables (or ' t ' ) these entries are 
ignored and a warning is printed. 



Example: 



t = numpy . linspace ( , 1, 500) 
x = numpy . cos ( . 2 *t ) 
y = numpy . sin ( . 2 *t ) 



histdic 




' t' : 


t 


' x' : 


X 


'y' : 


Y 



} 

dde . hist_f rom_arrays (histdic) 



hist_f rom_funcs (die, nn=101) 

Initialise the histories with the functions stored in the dictionary die. The keys are the variable names. The 
function will be called as f (t ) for t in [ -maxdelay , ] on nn samples in the interval. 

This function provides the simplest way to set the history. It is often convenient to use python lambda 
functions for f . This way you can define the history function in place. 

If any variable names are missing in the dictionaries, the history of these variables is set to zero and a 
warning is printed. If the dictionary contains keys not matching any variables these entries are ignored and 
a warning is printed. 

Example: Initialise the history of the variables x and y with cos and sin functions using a finer sampling 
resolution: 



from math import sin, cos 

histdic = { 

' x' : lambda t : cos ( . 2 *t ) , 
'y': lambda t: sin(0.2*t) 

} 

dde . hist_f rom_f uncs (histdic, 500) 

output_ccode ( ) 

run ( ) 

run the simulation 

class dde23 (eqns, params=None, noise=None, supportcode=", debug=False) 

This class translates a DDE to C and solves it using the Bogacki-Shampine method. 

Attributes of class instances: 

For user relevant attributes: 

self. sol Dictionary storing the solution (when the simulation has finished). The keys are the variable names 
and ' t ' corresponding to the sampled times. 

self.discont List of discontinuity times. This is generated from the occurring delays by propagating the 
discontinuity at t = 0. The solver will step on these discontinuities. If you want the solver to step onto 
certain times they can be inserted here. 

self.rseed Can be set to initialise the random number generator with a specific seed. If not set it is initialised 
with the time. 

self. hist Dictionary with the history. Don't manipulate the history arrays directly! Use the provided func- 
tions to set the history. 

self.Vhist Dictionary with the time derivatives of the history. 

For user less relevant attributes: 



selfdelays List of the delays occurring in the equations. 



self. chunk When arrays become to small they are grown by this number. 

self .spline Jck Dictionary which stores the tck spline representation of the solutions. (see 
scipy . interpolate) 

self.eqns Stores the eqn dictionary. 

selfparams Stores the parameter dictionary. 

selfsimul Dictionary of the simulation parameters. 

self. noise Stores the noise dictionary. 

self. debug Stores the debug flag. 

selfdelayhashs List of hashs for each delay (this is used in the generated C-code). 
selfvars List of variables extracted from the eqn dictionary keys. 
self. types Dictionary of C-type names of each variable. 
selfnptypes Dictionary of numpy-type names of each variable. 
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